# SIM CF version

set.seed(20190901)

plot_explicit_sd_units <- TRUE

iat_data_dummy <-
  iat_data %>%
  mutate(other = case_when(ethnicity == "other" ~ 1,
                           TRUE ~ 0),
         Ukr_and_Rus  = case_when(ethnicity == "Ukr and Rus" ~ 1,
                                  TRUE ~ 0),
         Russ_only= case_when(ethnicity == "Russian only" ~ 1,
                              TRUE ~ 0),
         Kherson = case_when(city == "Kherson" ~ 1,
                             TRUE ~ 0),
         Odesa  = case_when(city == "Odesa" ~ 1,
                            TRUE ~ 0),
         Kharkiv= case_when(city == "Kharkiv" ~ 1,
                            TRUE ~ 0),
         proRuss_vote = case_when(parl2014_3 == "Pro-Russia Voters" ~ 1,
                                  TRUE ~ 0),
         abstain = case_when(parl2014_3 == "Abstainers" ~ 1,
                             TRUE ~ 0))

##
## MODEL 9 (implicit bias as outcome variable)
##

m9_formula <- IAT ~ unemployed + family_in_Russia + russian_speaking_d + female +
  age_sd + proRuss_vote + abstain +
  other + Ukr_and_Rus + Russ_only + Kherson + Odesa + Kharkiv +
  other * Kherson +
  other * Odesa +
  other * Kharkiv +
  Ukr_and_Rus * Kherson +
  Ukr_and_Rus * Odesa +
  Ukr_and_Rus * Kharkiv +
  Russ_only * Kherson +
  Russ_only * Odesa +
  Russ_only * Kharkiv

m9_formula <- update.formula(old = m9_formula,
                             new = str_c(as.character(m9_formula)[2],
                                         "~",
                                         as.character(m9_formula)[3],
                                         "- proRuss_vote - abstain"))


# extract data
m9_data <- extractdata(m9_formula, iat_data_dummy)

# fit model, extract PEs and robust VCOV matrix
m9_dum_fit <- lm(m9_formula, data = m9_data)
m9_PEs <- m9_dum_fit$coefficients
m9_hcvoc <- vcovHC(m9_dum_fit)

# simulate coefficients
sims <- 10000
m9_simbetas <- mvrnorm(sims, m9_PEs, m9_hcvoc)

## MODEL 9 PREDICTED VALUES

# create counterfactual scenarios for expected values
scens <- 8
m9_scens <- cfMake(formula = m9_formula, data = m9_data, nscen = scens)

# change values of counterfactual scenarios
for(i in 1:scens) {
  m9_scens <- cfChange(m9_scens, "Kherson", x = 0, scen = i)
  m9_scens <- cfChange(m9_scens, "Odesa", x = 0, scen = i)
  m9_scens <- cfChange(m9_scens, "Kharkiv", x = 0, scen = i)
  m9_scens <- cfChange(m9_scens, "other", x = 0, scen = i)
  m9_scens <- cfChange(m9_scens, "Ukr_and_Rus", x = 0, scen = i)
  m9_scens <- cfChange(m9_scens, "Russ_only", x = 0, scen = i)
  m9_scens <- cfChange(m9_scens, "abstain", x = 0, scen = i)
  m9_scens <- cfChange(m9_scens, "proRuss_vote", x = 0, scen = i)
}

cities <- c(rep("Kharkiv", 2),
            rep("Kherson", 2),
            rep("Odesa", 2),
            rep("Kharkiv", 2))

for(i in 3:scens) {
  m9_scens <- cfChange(m9_scens, cities[i], x = 1, scen = i)
  
}

for(i in seq(2, 8, 2)) {
  m9_scens <- cfChange(m9_scens, "Russ_only", x = 1, scen = i)
  
}

# generate predicted values for counterfactual scenarios
m9_preds <- bind_cols(linearsimev(m9_scens, m9_simbetas, ci = 0.95))

# wrangle data to work w ggplot
colnames(m9_preds) <- c("fit", "lwr", "upr")

m9_data_ev_plot <-
  bind_cols(m9_scens$x,
            m9_preds) %>%
  mutate(ethnicity = case_when(Russ_only == 1 ~ "Russian only",
                               TRUE ~ "Ukrainian only"),
         city = case_when(Kherson == 1 ~ "Kherson",
                          Odesa == 1 ~ "Odesa",
                          Kharkiv == 1 ~ "Kharkiv",
                          TRUE ~ "Kyiv"),
         x_val = str_c(ethnicity, " X ", city),
         bias = "implicit bias") %>%
  dplyr::select(-proRuss_vote, -abstain, -other, -Ukr_and_Rus, -Russ_only,
                -Kherson, -Odesa, -Kharkiv)

## MODEL 9 FIRST DIFFERENCES

# create counterfactual scenarios for first differences
fd_scens <- 4
m9_fd_scens <- cfMake(formula = m9_formula, data = m9_data, nscen = fd_scens)

# change values of counterfactual scenarios
for(i in 1:fd_scens) {
  m9_fd_scens <- cfChange(m9_fd_scens, "Kherson", x = 0, xpre = 0, scen = i)
  m9_fd_scens <- cfChange(m9_fd_scens, "Odesa", x = 0, xpre = 0, scen = i)
  m9_fd_scens <- cfChange(m9_fd_scens, "Kharkiv", x = 0, xpre = 0, scen = i)
  m9_fd_scens <- cfChange(m9_fd_scens, "other", x = 0, xpre = 0, scen = i)
  m9_fd_scens <- cfChange(m9_fd_scens, "Ukr_and_Rus", x = 0, xpre = 0, scen = i)
  m9_fd_scens <- cfChange(m9_fd_scens, "Russ_only", x = 0, xpre = 0, scen = i)
  m9_fd_scens <- cfChange(m9_fd_scens, "abstain", x = 0, xpre = 0, scen = i)
  m9_fd_scens <- cfChange(m9_fd_scens, "proRuss_vote", x = 0, xpre = 0, scen = i)
}

cities <- c("Kyiv",
            "Kherson",
            "Odesa",
            "Kharkiv")

for(i in 2:fd_scens) {
  m9_fd_scens <- cfChange(m9_fd_scens, cities[i], x = 1, xpre = 1, scen = i)
  
}

for(i in 1:fd_scens) {
  m9_fd_scens <- cfChange(m9_fd_scens, "Russ_only", x = 1, scen = i)
  
}

# generate first differences for counterfactual scenarios
m9_fd_preds <- bind_cols(linearsimfd(m9_fd_scens, m9_simbetas, ci = 0.95))

# wrangle data to work w ggplot
colnames(m9_fd_preds) <- c("fit", "lwr", "upr")

m9_data_fd_plot <-
  bind_cols(m9_fd_scens$x,
            m9_fd_preds) %>%
  mutate(ethnicity = case_when(Russ_only == 1 ~ "Russian only",
                               TRUE ~ "Ukrainian only"),
         city = case_when(Kherson == 1 ~ "Kherson",
                          Odesa == 1 ~ "Odesa",
                          Kharkiv == 1 ~ "Kharkiv",
                          TRUE ~ "Kyiv"),
         x_val = str_c(ethnicity, " X ", city),
         bias = "implicit bias") %>%
  dplyr::select(-proRuss_vote, -abstain, -other, -Ukr_and_Rus, -Russ_only,
                -Kherson, -Odesa, -Kharkiv)

##
## MODEL 10 (EXPLICIT bias as outcome variable)
##

explicit_sd <- sd(iat_data$explicit)


m10_formula <- explicit ~ unemployed + family_in_Russia + russian_speaking_d + female +
  age_sd + proRuss_vote + abstain +
  other + Ukr_and_Rus + Russ_only + Kherson + Odesa + Kharkiv +
  other * Kherson +
  other * Odesa +
  other * Kharkiv +
  Ukr_and_Rus * Kherson +
  Ukr_and_Rus * Odesa +
  Ukr_and_Rus * Kharkiv +
  Russ_only * Kherson +
  Russ_only * Odesa +
  Russ_only * Kharkiv +
  IAT

# extract data
m10_data <- extractdata(m10_formula, iat_data_dummy)

# fit model, extract PEs and robust VCOV matrix
m10_dum_fit <- lm(m10_formula, data = m10_data)
m10_PEs <- m10_dum_fit$coefficients
m10_hcvoc <- vcovHC(m10_dum_fit)

# simulate coefficients
sims <- 10000
m10_simbetas <- mvrnorm(sims, m10_PEs, m10_hcvoc)

## MODEL 10 PREDICTED VALUES

# create counterfactual scenarios for expected values
scens <- 8
m10_scens <- cfMake(formula = m10_formula, data = m10_data, nscen = scens)

# change values of counterfactual scenarios
for(i in 1:scens) {
  m10_scens <- cfChange(m10_scens, "Kherson", x = 0, scen = i)
  m10_scens <- cfChange(m10_scens, "Odesa", x = 0, scen = i)
  m10_scens <- cfChange(m10_scens, "Kharkiv", x = 0, scen = i)
  m10_scens <- cfChange(m10_scens, "other", x = 0, scen = i)
  m10_scens <- cfChange(m10_scens, "Ukr_and_Rus", x = 0, scen = i)
  m10_scens <- cfChange(m10_scens, "Russ_only", x = 0, scen = i)
  m10_scens <- cfChange(m10_scens, "abstain", x = 0, scen = i)
  m10_scens <- cfChange(m10_scens, "proRuss_vote", x = 0, scen = i)
}

cities <- c(rep("Kharkiv", 2),
            rep("Kherson", 2),
            rep("Odesa", 2),
            rep("Kharkiv", 2))

for(i in 3:scens) {
  m10_scens <- cfChange(m10_scens, cities[i], x = 1, scen = i)
  
}

for(i in seq(2, 8, 2)) {
  m10_scens <- cfChange(m10_scens, "Russ_only", x = 1, scen = i)
  
}

# generate predicted values for counterfactual scenarios
m10_preds <- bind_cols(linearsimev(m10_scens, m10_simbetas, ci = 0.95))

# wrangle data to work w ggplot
colnames(m10_preds) <- c("fit", "lwr", "upr")
m10_data_ev_plot <-
  bind_cols(m10_scens$x,
            m10_preds) %>%
  mutate(ethnicity = case_when(Russ_only == 1 ~ "Russian only",
                               TRUE ~ "Ukrainian only"),
         city = case_when(Kherson == 1 ~ "Kherson",
                          Odesa == 1 ~ "Odesa",
                          Kharkiv == 1 ~ "Kharkiv",
                          TRUE ~ "Kyiv"),
         x_val = str_c(ethnicity, " X ", city),
         bias = "explicit bias") %>%
  dplyr::select(-proRuss_vote, -abstain, -other, -Ukr_and_Rus, -Russ_only,
                -Kherson, -Odesa, -Kharkiv)


## MODEL 10 FIRST DIFFERENCES

# create counterfactual scenarios for expected values
fd_scens <- 4
m10_fd_scens <- cfMake(formula = m10_formula, data = m10_data, nscen = fd_scens)

# change values of counterfactual scenarios
for(i in 1:fd_scens) {
  m10_fd_scens <- cfChange(m10_fd_scens, "Kherson", x = 0, xpre = 0, scen = i)
  m10_fd_scens <- cfChange(m10_fd_scens, "Odesa", x = 0, xpre = 0, scen = i)
  m10_fd_scens <- cfChange(m10_fd_scens, "Kharkiv", x = 0, xpre = 0, scen = i)
  m10_fd_scens <- cfChange(m10_fd_scens, "other", x = 0, xpre = 0, scen = i)
  m10_fd_scens <- cfChange(m10_fd_scens, "Ukr_and_Rus", x = 0, xpre = 0, scen = i)
  m10_fd_scens <- cfChange(m10_fd_scens, "Russ_only", x = 0, xpre = 0, scen = i)
  m10_fd_scens <- cfChange(m10_fd_scens, "abstain", x = 0, xpre = 0, scen = i)
  m10_fd_scens <- cfChange(m10_fd_scens, "proRuss_vote", x = 0, xpre = 0, scen = i)
}

cities <- c("Kyiv",
            "Kherson",
            "Odesa",
            "Kharkiv")

for(i in 2:fd_scens) {
  m10_fd_scens <- cfChange(m10_fd_scens, cities[i], x = 1, xpre = 1, scen = i)
  
}

for(i in 1:fd_scens) {
  m10_fd_scens <- cfChange(m10_fd_scens, "Russ_only", x = 1, scen = i)
  
}

# generate predicted values for counterfactual scenarios
m10_fd_preds <- bind_cols(linearsimfd(m10_fd_scens, m10_simbetas, ci = 0.95))

# wrangle data to work w ggplot
colnames(m10_fd_preds) <- c("fit", "lwr", "upr")

m10_data_fd_plot <-
  bind_cols(m10_fd_scens$x,
            m10_fd_preds) %>%
  mutate(ethnicity = case_when(Russ_only == 1 ~ "Russian only",
                               TRUE ~ "Ukrainian only"),
         city = case_when(Kherson == 1 ~ "Kherson",
                          Odesa == 1 ~ "Odesa",
                          Kharkiv == 1 ~ "Kharkiv",
                          TRUE ~ "Kyiv"),
         x_val = str_c(ethnicity, " X ", city),
         bias = "explicit bias") %>%
  dplyr::select(-proRuss_vote, -abstain, -other, -Ukr_and_Rus, -Russ_only,
                -Kherson, -Odesa, -Kharkiv)

##
## COMBINE BOTH PLOTS FOR EXPECTED VALUES
##


m9_m10_ev <-
  m9_data_ev_plot %>%
  bind_rows(m10_data_ev_plot)

m9_m10_ev$city <-
  factor(m9_m10_ev$city,
         levels = c("Kyiv",
                    "Kherson",
                    "Odesa",
                    "Kharkiv"))


##
## COMBINE M9 AND M10 FIRST DIFF PLOTS
##

m9_m10_fd <-
  m9_data_fd_plot %>%
  bind_rows(m10_data_fd_plot)

m9_m10_fd$city <-
  factor(m9_m10_fd$city,
         levels = c("Kyiv",
                    "Kherson",
                    "Odesa",
                    "Kharkiv"))


m9_m10_fd$ethnicity <-
  factor(m9_m10_fd$ethnicity,
         levels = c("Ukrainian only",
                    "other",
                    "Ukr and Rus",
                    "Russian only"))


##
##  M9 AND M10 EXPECTED VALUES FOR ALL CITIES, ETHNICITIES, PARL VOTE CATS
##

# create counterfactual scenarios for expected values
all_scens <- 48

other <- c(rep(1, 12), rep(0, 36))
ukr_rus <- c(rep(0, 12), rep(1, 12), rep(0, 24))
rus_only <- c(rep(0, 24), rep(1, 12), rep(0, 12))

kherson <- c(rep(c(rep(1, 3), rep(0, 9)), 4))
odesa <- c(rep(c(rep(0, 3), rep(1, 3), rep(0, 6)), 4))
kharkiv <- c(rep(c(rep(0, 6), rep(1, 3), rep(0, 3)), 4))

abstain <- c(rep(c(1, 0, 0), 16))
pro_russ <- c(rep(c(0, 1, 0), 16))

m9_scens_all <- cfMake(formula = m9_formula, data = m9_data, nscen = all_scens)

# change values of counterfactual scenarios
for(i in 1:all_scens) {
  m9_scens_all <- cfChange(m9_scens_all, "Kherson", x = kherson[i], scen = i)
  m9_scens_all <- cfChange(m9_scens_all, "Odesa", x = odesa[i], scen = i)
  m9_scens_all <- cfChange(m9_scens_all, "Kharkiv", x = kharkiv[i], scen = i)
  m9_scens_all <- cfChange(m9_scens_all, "other", x = other[i], scen = i)
  m9_scens_all <- cfChange(m9_scens_all, "Ukr_and_Rus", x = ukr_rus[i], scen = i)
  m9_scens_all <- cfChange(m9_scens_all, "Russ_only", x = rus_only[i], scen = i)
  m9_scens_all <- cfChange(m9_scens_all, "abstain", x = abstain[i], scen = i)
  m9_scens_all <- cfChange(m9_scens_all, "proRuss_vote", x = pro_russ[i], scen = i)
}

# generate predicted values for counterfactual scenarios
m9_preds_all <- bind_cols(linearsimev(m9_scens_all, m9_simbetas, ci = 0.95))

# wrangle data to work w ggplot
colnames(m9_preds_all) <- c("fit", "lwr", "upr")
m9_all_data_ev_plot <-
  bind_cols(m9_scens_all$x,
            m9_preds_all) %>%
  mutate(ethnicity = case_when(Russ_only == 1 ~ "Russian only",
                               other == 1 ~ "other",
                               Ukr_and_Rus == 1 ~ "Ukr and Rus",
                               TRUE ~ "Ukrainian only"),
         city = case_when(Kherson == 1 ~ "Kherson",
                          Odesa == 1 ~ "Odesa",
                          Kharkiv == 1 ~ "Kharkiv",
                          TRUE ~ "Kyiv"),
         parl2014_3 = case_when(proRuss_vote == 1 ~ "Pro-Russia Voters",
                                abstain == 1 ~ "Abstainers",
                                TRUE ~ "Anti-Russia Voters"),
         x_val = str_c(ethnicity, " X ", city),
         bias = "implicit bias") %>%
  dplyr::select(-proRuss_vote, -abstain, -other, -Ukr_and_Rus, -Russ_only,
                -Kherson, -Odesa, -Kharkiv)

m10_scens_all <- cfMake(formula = m10_formula, data = m10_data, nscen = all_scens)

# change values of counterfactual scenarios
for(i in 1:all_scens) {
  m10_scens_all <- cfChange(m10_scens_all, "Kherson", x = kherson[i], scen = i)
  m10_scens_all <- cfChange(m10_scens_all, "Odesa", x = odesa[i], scen = i)
  m10_scens_all <- cfChange(m10_scens_all, "Kharkiv", x = kharkiv[i], scen = i)
  m10_scens_all <- cfChange(m10_scens_all, "other", x = other[i], scen = i)
  m10_scens_all <- cfChange(m10_scens_all, "Ukr_and_Rus", x = ukr_rus[i], scen = i)
  m10_scens_all <- cfChange(m10_scens_all, "Russ_only", x = rus_only[i], scen = i)
  m10_scens_all <- cfChange(m10_scens_all, "abstain", x = abstain[i], scen = i)
  m10_scens_all <- cfChange(m10_scens_all, "proRuss_vote", x = pro_russ[i], scen = i)
}

# generate predicted values for counterfactual scenarios
m10_preds_all <- bind_cols(linearsimev(m10_scens_all, m10_simbetas, ci = 0.95))

if(plot_explicit_sd_units) {
  m10_preds_all <- m10_preds_all / explicit_sd
}

# wrangle data to work w ggplot
colnames(m10_preds_all) <- c("fit", "lwr", "upr")
m10_all_data_ev_plot <-
  bind_cols(m10_scens_all$x,
            m10_preds_all) %>%
  mutate(ethnicity = case_when(Russ_only == 1 ~ "Russian only",
                               other == 1 ~ "other",
                               Ukr_and_Rus == 1 ~ "Ukr and Rus",
                               TRUE ~ "Ukrainian only"),
         city = case_when(Kherson == 1 ~ "Kherson",
                          Odesa == 1 ~ "Odesa",
                          Kharkiv == 1 ~ "Kharkiv",
                          TRUE ~ "Kyiv"),
         parl2014_3 = case_when(proRuss_vote == 1 ~ "Pro-Russia Voters",
                                abstain == 1 ~ "Abstainers",
                                TRUE ~ "Anti-Russia Voters"),
         x_val = str_c(ethnicity, " X ", city),
         bias = "explicit bias") %>%
  dplyr::select(-proRuss_vote, -abstain, -other, -Ukr_and_Rus, -Russ_only,
                -Kherson, -Odesa, -Kharkiv)

m9_m10_ev_all <-
  m9_all_data_ev_plot %>%
  bind_rows(m10_all_data_ev_plot)

m9_m10_ev_all$city <-
  factor(m9_m10_ev_all$city,
         levels = c("Kyiv",
                    "Kherson",
                    "Odesa",
                    "Kharkiv"))

m9_m10_ev_all$parl2014_3 <-
  factor(m9_m10_ev_all$parl2014_3,
         levels = c("Anti-Russia Voters",
                    "Pro-Russia Voters",
                    "Abstainers"))

m9_m10_ev_all$ethnicity <-
  factor(m9_m10_ev_all$ethnicity,
         levels = c("Russian only",
                    "Ukr and Rus",
                    "other",
                    "Ukrainian only"))

ev_full_plot <-
  ggplot(data = m9_m10_ev_all) +
  geom_pointrange(aes(x = ethnicity, y = fit, ymin = lwr, ymax = upr,
                      # linetype = bias,
                      color = bias),
                  position = position_dodge(width = -.5)) +
  geom_hline(yintercept = 0, color = "gray", linetype = "dashed", size = .2) +
  scale_color_manual(values = c("dark gray", "black")) +
  # scale_y_continuous(breaks = seq(-1.5, 1.5, .5),
  #                    limits = c(-1.7, 1.7)) +
  facet_grid(parl2014_3 ~ city) +
  labs(#title = "",
    x = NULL,
    y = "predicted bias (values greater than zero are pro-Ukraine)") +
  theme_tufte() +
  theme(legend.position = "bottom",
        legend.title = element_blank()) +
  coord_flip()

## PREDICTED PROBS AND FIRST DIFFS FOR RUSSIAN ONLY

# run_sim_cf(plot_explicit_sd_units = TRUE)

# generate Figure 4; generate Figure E.4

exp_ev <-
  ggplot(data = m9_m10_ev %>%
           filter(ethnicity == "Russian only" &
                    bias == "explicit bias")) +
  geom_pointrange(aes(x = city, y = fit, ymin = lwr, ymax = upr),
                  color = "dark gray",
                  position = position_dodge(width = -.5)) +
  geom_hline(yintercept = 0,  linetype = "dashed", size = .2) +
  # scale_color_manual(values = c("dark gray", "black")) +
  scale_y_continuous(#position = "right", 
    limits = c(-1, 3.5)) +
  # facet_wrap( ~ bias, nrow = 1) +
  labs(title = "Explicit bias (ethnic Russians)",
       x = NULL,
       y = "predicted bias") +
  theme_tufte() +
  theme(legend.position = "none",
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5)) 

imp_ev <-
  ggplot(data = m9_m10_ev %>%
           filter(ethnicity == "Russian only" &
                    bias == "implicit bias")) +
  geom_pointrange(aes(x = city, y = fit, ymin = lwr, ymax = upr),
                  color = "black",
                  position = position_dodge(width = -.5)) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .2) +
  scale_color_manual(values = c("dark gray", "black")) +
  scale_y_continuous(position = "right", 
                     limits = c(-.22, .77)) +
  # facet_wrap( ~ bias, nrow = 1) +
  labs(title = "Implicit bias (ethnic Russians)",
       x = NULL,
       y = "predicted bias") +
  theme_tufte() +
  theme(legend.position = "none",
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5)) 

exp_fd <-
  ggplot(data = m9_m10_fd %>%
           filter(bias == "explicit bias")) +
  geom_pointrange(aes(x = city, y = fit, ymin = lwr, ymax = upr),
                  color = "dark gray",
                  position = position_dodge(width = -.5)) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .2) +
  scale_color_manual(values = c("dark gray", "black")) +
  scale_y_continuous(limits = c(-4, 0.8)) +
  # facet_wrap( ~ bias, nrow = 1) +
  labs(title = "Explicit bias change\n(first difference, Russ. vs Ukr.)",
       x = NULL,
       y = "change in bias") +
  theme_tufte() +
  theme(legend.position = "none",
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5))

imp_fd <-
  ggplot(data = m9_m10_fd %>%
           filter(bias == "implicit bias")) +
  geom_pointrange(aes(x = city, y = fit, ymin = lwr, ymax = upr),
                  color = "black",
                  position = position_dodge(width = -.5)) +
  geom_hline(yintercept = 0, linetype = "dashed", size = .2) +
  # scale_color_manual(values = c("dark gray", "black")) +
  scale_y_continuous(position = "right",
                     limits = c(-0.7, 0.14)) +
  # facet_wrap( ~ bias, nrow = 1) +
  labs(title = "Implicit bias change\n(first difference, Russ. vs Ukr.)",
       x = NULL,
       y = "change in bias") +
  theme_tufte() +
  theme(legend.position = "none",
        legend.title = element_blank(),
        plot.title = element_text(hjust = 0.5))

ev_fd_rus_plot_name <-"./plots/ev_fd_plot_n_no_parl.pdf"


ggsave(filename = ev_fd_rus_plot_name,
       arrangeGrob(exp_ev, imp_ev, exp_fd, imp_fd, ncol = 2),
       width = 7, height = 7)


ev_all_plot_name <- "./plots/ev_plot_all_no_parl.pdf"


# generate Fig F.5
ggsave(ev_full_plot,
       filename = ev_all_plot_name,
       width = 7, height = 7)






